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ABSTRACT 


Joint inversion of magnetic, Earth rotation, geoid, and seismic 
data for a unified model of the coupled core-mantle system is proposed 
and shown to be possible. A sample objective function is offered and 
simplified by targeting results from independent inversions and summary 
travel time residuals instead of original observations. These "data" 
are parameterized in terms of a very simple, closed model of the 
topographically coupled core-mantle system. Minimization of the 
simplified objective function leads to a non-linear inverse problem; an 
iterative method for solution is presented. Parameterization and method 
are emphasized; numerical results are not presented. 
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1. INTRODUCTION 

Geophysicists working with different types of data are probing 
Earth's deep interior (see, e.g., Lay, 1989). For example, geomagnetic 
data have been used to estimate fluid motions near the top of the core 
(Ball, Kahle & Vestine, 1969; Voorhies, 1984, 1986a, b, 1988; LeMouSl , 
Gire & Madden, 1985; Whaler & Clarke, 1988); seismic data have been used 
to estimate laterally heterogeneous mantle structure and core-mantle 
boundary - hereafter denoted CMB - topography (Morel li & Dziewonski, 
1987); gravity and geodetic data have been combined with seismic 
estimates of Earth structure to estimate CMB topography (Hager et al . , 
1985); and estimates of surficially geostrophic core motions have been 
combined with estimates of CMB topography to calculate the topographic 
torque exerted by the core on the mantle and the implied changes in 
"solid" Earth rotation (Speith et al., 1986). The latter uses results 
from independent or "disjoint" inversions of different geophysical data 
types to forwardly model decade fluctuations in solid Earth rotation. 

I propose joint inversion of diverse geophysical data types for a 
unified model of the coupled core-mantle system. The plan merges 
magnetic, Earth rotation, geoid, and seismic data into one objective 
function which, when suitably weighted, constrained, and parameterized, 
can be minimized with respect to the parameters of a unified deep Earth 
model. The goal is to develop, parameterize, and test hypotheses about 
Earth's deep interior against all relevant types of data. 

Curiously, the philosophical foundation for this type of inversion 
has been questioned. Clearly, much can be learned from experiments 
designed to isolate those data which are thought to be most sensitive to 
some particular property of the Earth. This approach can yield decisive 
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tests of particular hypotheses; yet one need not always lose sight of 
the forest for the trees. Some properties of Earth's deep interior 
(e.g., CMB topography) can contribute signals to many kinds of data 
yet are apparently not uniquely determined by any single kind of data. 

In such cases, more plausible estimates of the properties might be 
obtained by using more than one kind of data. 

To do so, a merged data set may be compiled and used to estimate 
parameters of models of the Earth properties. One can hypothesize that 
signals from properties which are not modeled, and from parameters which 
are not estimated, do not vastly exceed the residuals indicated by a 
weighted least-squares fit of the modeled parameters to the data. This 
hypothesis can, in turn, be investigated by fitting more data and more 
types of data with more complete models of more Earth properties. 

To this end, I offer a sample objective function and parameterize 
it in terms of a simplified, mechanically coupled, core-mantle system. 
The sample "data" considered are slowly varying geomagnetic potential 
coefficients, decade fluctuations in the angular velocity of the solid 
Earth, static gravitational potential coefficients, and summary seismic 
travel time residuals relative to a laterally homogeneous Earth model. 
The parameters describe a piecewise steady core surface velocity field, 
a perturbation density field in the mantle, and a CMB topography 
function. The system is closed by supposing surficially, indeed 
tangentially (Backus & LeMoufil , 1986), geostrophic core motions and 
relations between perturbation seismic wave speeds and perturbation 
density in the mantle. Even for this simple Earth model, minimization 
of the sample objective function leads to a non-linear inverse problem; 
an iterated, linearized method of solution is presented. 
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This sample is intended to provide a foundation for more realistic 
deep Earth models which might include: a superior mean state; mantle 
dynamics and rheology; richer core dynamics; magnetic, viscous, and 
gravitational core-mantle coupling; and thermal and compositional core- 
mantle interactions. More work will be needed on the problems of how to 
parameterize such models, include more kinds of data (e.g., free 
oscillations and plate motions), and apply more constraints (e.g., from 
mineral physics and low-frequency gravity and deformation studies); and 
on problems of uniqueness, accuracy, and method. 

2. AN OBJECTIVE FUNCTION SIMPLIFIED 
Let r be the position vector in geobarycentric spherical polar 
coordinates radius r, colatitude 0 , and east longitude let t be time; 
and let observational data and Earth model predictions be denoted 
respectively by d and p subscripts on the following variables: 

B is the geomagnetic flux density vector; 

Q is the apparent angular velocity vector at the surface of the 
solid Earth, technically including plate motions; 
g is the gravitational acceleration vector; 

T is the travel time of seismic phase £ from the source at ( r ' , t ' ) 
to the receiver at (r,t); 

W* is the weight function which is taken to be the inverse squared 
uncertainty in datum of type x at (r,t), but which can be 
generalized to a weight matrix for discrete data on the 
expectation of correlated errors; and 
Xi are multipliers giving weight and proper units to the 
Ki which represent geophysical constraints (e.g., small density 
perturbations, finite surficially geostrophic core fluid 
velocity (Mach number < 1), smooth CMB topography, etc.). 
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One suitable objective function to be minimized in a joint 
inversion of geophysical data is 

A? = {the weighted residual variance relative to the 

magnetic + Earth rotation + geoid + seismic data} 

+ {other geophysical constraints} 

= Am^ + Aer^ + Age^ + As^ + XjAKi^ 

or, in the foregoing notation, 

r 0 2r r tf 0 m 

A 2 = J / J J { [Bd(r,t) - B p (r,t)]2w (r,t) 
ri 0 0 t 0 . . „ er 

+ [Qd(r,t) - D p (r,t)]2w (r,t) 

+ [gd (r, t) - g p (r,t)]2w 9e (r,t) 

r'o 2 jt jt t'f 

+ [/ Iff [Td(^;r' ,t';r,t) - Tp(^;r’ ,t' ;r,t)]2 
r'-jOOt'os 

W (f;r' ,t' ,r,t)r'2sin0'dt'd0'd^'dr'] 

+ [XiKi (r,t)2] }r2sin0dtd0d^dr . 

It is convenient to view the results of measurements over small portions 

of the 4-volume of integration as discrete data. Then 

„ merges 

A2 = [dj - p j ] Wjk [dk - pk] + <X i K^> 

where dj is an element of the merged data vector, pj is an element of 
the merged prediction vector, and repeated subscripts denote summation. 

This objective function could be generalized to include other types 
of data; yet it already seems too ambitious and the data have already 
been reduced and analyzed "disjointly" . The suggestion is to build upon 
this rich tradition by replacing the diverse data types with either more 
tractable models thereof or data residuals relative to such models. One 
such approach begins with the following, de-subscripted, "data": 

(1) slowly varying, broad-scale spherical harmonic models B of the 
observed geomagnetic field (e.g., IAGA 1988); 
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(2) length-of-day and polar motion data corrected as possible for 
nutation, precession, and tidal effects, either low-pass 
filtered or corrected for atmospheric and hydro-cryospheric 
effects (e . g . , Stephenson & Morrison, 1984), and then fitted 
with, say, a piecewise linear function D(t); 

(3) broad-scale spherical harmonic models of the steady part of the 
gravitational field g (e.g., Marsh et al . , 1988) - preferably 
corrected for surface topographic and crustal sources; and 

(4) summary travel time residuals T relative to a laterally 
homogeneous seismic Earth model (e.g., Dziewonski & Anderson, 
1981) which specifies the axisymmetric mean state 

(Vpo. VSo,/?o.PorK 0 ,/*o) on reference ellipsoids (or spheres). 
Effects of external fields on (1) and (3) are small. Effects of plate 
motions on (2) are omitted for now; the piecewise linear Q(t) fitted to 
corrected, low-pass filtered Od should capture the decade fluctuations 
of interest here. Summary travel time residuals in (4) are averaged 
over closely spaced ray paths (Creager & Jordan, 1986) to reduce effects 
of small-scale structure, oversampling, and colinearity. Alternately, 
one might use (4b) a model of laterally heterogeneous phase speeds Vp 
and Vs, or (4c) spherical harmonic models of the travel time residuals 
at all (summary) receivers for each phase from each (summary) source. 

Let the reference surface of mean radius lal = a 2 6.3712 Mm 
enclose the internal sources of scaloidal B and g. On and above this 
surface we have 

Bd = B(r, t) + <5b(r,t) B = -VV 

with internal scalar magnetic potential 
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w a n+1 n r in, % m _m 

V = a E [-] E [g n (t)cosm^ + h n (t)sinmf]Pn(cos0) 
n=l r m=0 

and radial magnetic component 
9V 

B r (a,t) = = E gi(a r t) Si(0,^) = gi Si = g T S . (1) 

9r i " ' 

Here the (gn m .hn m ) are the Gauss coefficients, P n m is the Schmidt- 

normalized associated Legendre polynomial of degree n and order m, gi is 

an element of the ordered column vector g of radial magnetic field 

coefficients, Si is an element of the ordered vector S of spherical 

harmonics of degree n(i) and order m(i), and a T superscript indicates 

the transpose (Voorhies, 1986b). Moreover, 

Od(t) = Q 0 + 0(t) + $i#(t) (2) 

where Q 0 is the mean angular velocity of the solid Earth and Q(t) is the 
piecewise linear representation of the decade fluctuations. Furthermore, 
with g 0 (r) = -VU 0 being the mean gravity caused by the mean density 
p 0 {r) , and with tidal and other time-varying effects represented by <5g, 

9d = 9o(r) + g(r) + <5g(r,t) . 

Steady g(r) = -VU has steady perturbation gravitational potential 

GM E w a n+1 n m m m 

U = — E [-] E [c n cosm^ + s n sinm^]P n (cos0) 
a n=0 r m=0 

and steady perturbation radial gravitational component 
9U 

gr (a) = = E ci (a) Siie.f) = ciSi = cTs . (3) 

9r i 

Here G is the gravitational constant, Me is Earth's mass, (c n m ,s n m ) are 
the steady perturbation gravitational potential coefficients, and ci is 
an element of the ordered vector c of steady perturbation radial 
gravitational field coefficients. Finally, 
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T(<r»r' ,t' ;r,t) = Td(f;r' ,t' ;r,t) - T 0 (£;r' ,t' ;r,t) = Ti = Tdi - T 0 i (4) 
where, for each seismic phase $ and each (summary) source-receiver pair 
denoted by i = i (£;r' ,t' ;r,t) , the (summary) travel time residual Ti is 
relative to the travel time predicted by the reference Earth model T 0 i . 
Hopefully gi, fl i , c i , and Ti would be corrected for external and crustal 
effects; such corrections to low-degree geomagnetic main field models, 
piecewise linear fits to low-pass filtered Earth rotation data, a low- 
degree gravity model, and perhaps summary travel times should be small. 

With the foregoing reduction of data types, and targeting the 
evolution of the radial magnetic component and the steady perturbation 
radial gravitational component on a, the redefined objective function is 

„ tf 2ff m 

= / j j [B r (a , t ) - B rp (a,t)]2w (a,t)sin0d0d^dt 
t 0 0 0 

tf er 

+ J C0(t) - 0 p (t)]2w (t)dt 

to 

2ir jt ge 

+ J J [gr(a) - grp(a)] 2 W (a)sin0d0d^ 

0 0 

+ [Ti - T pi ]W S ij[Tj - T pj ] + <XjKi(r,t)2> . (5) 

Scalar weight functions W™ and W9 e can be derived from the (increasingly 
realistic) error covariance matrices for the geopotential field models, 
while W^r should reflect uncertainty and error in the piecewise linear 
fit Q(t). The diagonal elements of the matrix W s should reflect the 
standard deviation of each summary travel time residual; hypocenter 
uncertainties may suggest non-trivial off-diagonal elements. With eqns. 
(1) and (3), B r p(a,t) = 7 i S i , g r p(a) = fiSi, the three components of 
Op ( t ) written tl i p ( t ) = u\ (t) , and with T p |< = rk, evaluation of the 
weighted surface integrals in eq. (5) yields 
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tf m 

A 2 = / C(9i - 7i) w i j (9j " 7j) Jdt 

to 

tf er 

+ J [ (Q i - wi) Wij (flj - wj)]dt 
to 

+ [(ci - £i) w ij ( c j " 

+ [(Ti - ri) Wij (Tj - rj)] + <XiKi(r f t)2> . (6) 

Other types of geophysical data may be included by adding suitably 
weighted terms to the objective function given by eq. (6). 

Select changes in the weight matrices of eq. (6) can be made so as 
to transform the objective function. For example: (i) geomagnetic and 
gravity coefficients can be fitted through degree and order 10 without 
biasing the higher degree 7i and £i by setting WW-jj and W9 e ij equal to 
zero for either i or j greater than 120; (ii) replacing W^ij with the 
inverse of the error covariance matrix for the radial magnetic 
coefficients targets the scalar geomagnetic potential instead of the 
radial component alone; or (iii) W 5 -jj might be culled to restrict 
attention to particular phases like PKP or ScS. 

3. PARAMETERIZATION 

The parameterization of predictions 71, wj, fk. and r 1 offered here 
is based on a very simple Earth model with the following attributes. 

(1) Geomagnetic secular variation is attributed to piecewise steady 
(Voorhies & Backus, 1985) frozen-flux motional induction by a 
tangentially geostrophic (Backus & LeMou&l , 1986) fluid velocity 
field v s at the top of a roughly spherical core of mean radius let 
= c £ 3.48 Mm (no subscript or underscore) which is surrounded by a 
comparatively rigid and magnetically source-free mantle. 
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(2) Decade fluctuations in the angular velocity of the solid Earth are 

attributed to the mechanical torque L exerted by the perturbation 
pressure p'(c,0,j»;t) associated with v s on the topography h(0,/) of 
the CMB. The CMB is the locus of points r c : r c = c + h(0,^). 

(3) The perturbation gravity field is attributed to the perturbation 

density in the solid Earth p' (r£c+h,0,^) and the effect of CMB 
topography h. Density perturbations within the core are omitted. 

(4) Seismic travel time residuals are due to h, p' , and perturbations in 

the bulk and shear moduli K and p. 

3.1 Magnetic 

In this very simple model, within each subinterval during which 
steady flow is presumed the time rate of change of the predicted radial 
field at the CMB is, to order zero in h, 

Brp( c «t) = Vs* [Brp(c.t) v s (c)] = TjSj 

where the over-dot indicates the partial time derivative and Tj is an 
element of the ordered vector of predicted radial magnetic field 
coefficients at c. With the diagonal upward continuation matrix Tij 
having elements (c/a) n (i) + 2 for harmonic degree n(i), 

Brp( a <t) ■ 7 i S i = (TijTjjSi . 

The surficial fluid velocity is expressed in terms of the streamfunction 
-T m and the velocity potential -U m 

v s (c) = rxV s T m + V s U m Tm = fli $i U"i = p\S\ . 

These expressions imply 

7m = Tmk{[XiXi jk]®j + [ri Y ijk]/ , j} 

= Tmk{Xi jk^ioj + Yi j kf i/?j } = A m lv] (7) 
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where, in the last step, the sums over i and k have been performed and 
is an element of the concatenated vector of streamfunction 
coefficients aj and velocity potential coefficients Pj. Equation (7) is 
but 7 = Av. The time-varying elements A m i of matrix A depend on r i , 
hence Bpp, and thus upon the velocity field coefficients v]. The 
inverse problem is therefore non-linear; for the iterative linearized 
approach suggested in section 4, the A m l are first calculated from the 
gi and are recalculated on each deep iteration (Voorhies, 1987a, b, 

1988). Formulae for Xijk = -Xjik, Yijk. and A m l are given elsewhere for 
the linear case (Voorhies, 1986b); formulae for the non-linear case have 
been presented (Voorhies, 1987a), posted (Voorhies, 1988), and are in 
typescripts (available by request) detailing the steady core flow 
estimation methods applied routinely at GSFC. 

Tangential geostrophy (Ball et al . , 1969; Backus & LeMoubl , 1986) 
seems awkward to enforce. In contrast, it is easy to damp departures 
from a geostrophic radial vorticity balance (Voorhies, 1986b, c). Then 
V s *[v s cos0] 2 0 and downwelling implies poleward flow (Voorhies, 1987c). 
Such flows are but "surficially geostrophic" (Voorhies, 1990); 
subsequent supposition of tangential geostrophy allows calculation of 
the perturbation pressure field on the sphere c from the vi . Steady 
perturbation pressure "maps" so derived at GSFC show fair agreement with 
those derived from the work of C. Gire and J.-L. LeMou&l - the first, I 
believe, to produce such maps (D. Jault, 1989, personal communication). 
3.2 Earth Rotation 

The reference mantle has steady principle moments of inertia 
(A, B-A, C) in geobarycentric Cartesian coordinates (x,y,z) with Q 0 
parallel to the z axis; the time rate of change of the predicted angular 


10 



velocity vector is, according to the Euler equations, 

Wx = [Lx “ (Ho + w z ) (C-B)wy]/A - [Lx - Qo(C - A)wy]/A 

Wy = [Ly + (n 0 + w z ) (C-A)w x ]/A z [Ly + Qo(C-A)«x]/A 

w z = [L z - «x w y (B"A)]/C — L z /C 

with the approximations good to first order in I vi I /0 o « 1. Voorhies 
(1990) shows that, to first-order accuracy in the asphericity Ihl/c « 1 
and in I V s h I « 1, the topographic torque is (omitting Lorentz terms) 

2jt i 

L = 20 0 PcC 3 / / [h (0 , ^) v s ( r c =c ,6,f) cos0] s i n0d0d^ 

0 0 

With h(0,^) = hjSi one has, to first-order accuracy, 

Lk = [hiQ*ijk]vj = Q*i jk h i v j 

«k(t) = Qi jk^i vj + Ekj«j(t) = Qi jk^i vj + qk(t) (8) 

• • 

or simply u = Q T hv + q. Because wk(t) depends on piecewise steady hv s , 

hence h i vj , and because qk(t) depends on tfj(t), the inverse problem is 

non-linear. In the iterative linearized approach (see section 4), the 

qk, and perhaps A, B*A, C and off diagonals of the inertia tensor, 

should be recalculated on each deep iteration. Then it is convenient to 

introduce the matrices Q^kj = Qijk h i and Q h ki = Qijkvj- 

With initial conditions gi(a,t 0 ) = 7 i(a,t 0 ) and fij(t 0 ) = wj(t 0 ), 

and simply supposing steady flow from to to tf, the magnetic and Earth 

rotation portion of eq. (6) is 

tf m er 

/ [ (9i - 7i) w ij(9j - 7j)l + t( n i ~ «i) w ij(°j ~ w j)] dt 

to 

tf t . . m t . 

- / {[/ (91 - 7 i)dt']Wij[J (gj - 7j)df] + 

to t 0 ti 
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t * . er t • • 

[/ (Oi - «i)dt']Wij[/ (flj - wj)dt']}dt 
to to 

or, by eqns. (7) and (8) , 


= [g 


T m • 

Av] W [g - Av] + [0 


T T er • T 

Q hv - q] W [fi - Q hv - q] 


( 9 ) 


where the number of underscores denotes tensor rank, the first over-bar 
indicates dummy time integration from to to t, and the second over-bar 
indicates time integration from t 0 to tf. Besides the explicit non- 
linear dependence on hv, there are non-linearities implicit in 
and g(w(h, v) ) • 

3.3 Geoid 

The suggested geoid parameterization is in terms of spherical 
harmonic coefficients for a perturbation density which is independent of 
radius within each of K layers 


p' (rk<r*r|<+l,0,^) = |>(rk<r£rk+l)] i Si [9,f) = ^kiSi 

and the mean state density contrast across the CMB, bp = pc - p m, times 
the CMB topography 


p' (fo-c ,9 ,#) ~ A/?hiSi . 

The perturbation geopotential is 
P'{r') n 

U(r) = -GfJJ r^sinfl'dA'd^'dr' . 

Ir - r' I 

Now g r = -U, r ; with lal = Irl > Ir'l, expanding I r-r ’ I -1 yields 

a 2ir n r ' n(j)+2 

g r (a) = -G / / / {[/>' (r' )] iSi (*' )}{[n(j)+l] (-) 

c 0 0 a 

Sj {6’ ,(j)' )Sj (fl.^Jsinfl'dfl'd^'dr 1 
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a 4jG[n(j)+l] r > n(j)+2 

= -/ [p' (' r ')]i{ — " — — ( ) 5ij}Sj(0 ( ^)dr' 

c [2n(j)+l] a 

K -4*-aG [n(j)+l] rk+l n(i)+3 r|< n(i)+3 

k=l [2n(j)+l] [n(i)+3] a a J J 

-4iraG [n(j)+l] c n(i)+3 

+ A/?hi{ (-) Sij}Sj 

[2n(j)+l] [n(i)+3] a 

= [/^kiGikjJSj + [hiHijjSj = [^kiGikj + H i j h i ] S j 
where <5ij is the Kroenecker delta. Because Gikj (and Hij) are zero for 
i*j, the connection between p' and gr(a) is harmonically pure (as is 
that between h and g r (a)). It is convenient to reorder the ^kj into a 
single vector p i with 1 = (k-j ) j max + j so that />kiGikj ■ Fj i/?i and 


9rp(«) - [Fjl/»1 + Hj i hi ] Sj = fjSj = £ T S . (10) 

The geoid contribution to eq. (6) is thus 

[(ci - fi) W?j (cj - fj)] 

ge 

= t c i " Fi i /7i - Himhni}Wi j{cj - Fji/j] - Hj m h m } 
ge 

= [c - F p - Hh]Tw [c - F p - Hh] . (II) 

3.4 Seismic 

The predicted (summary) travel time perturbation ri is the 
perturbation slowness ^ i ( r) integrated along the ray path L(i) for phase 
f(i) from source r'(i) to receiver r(i) 
r(i) i 

r\ = J * (r) d L i ' 
r' (i) 

No sum is performed over superscripts. For local (P or S) phase speed 
V = V 0 + A V = (s 0 + ji) - l, s 0 - 1/V 0 and y z -AV/V 0 2 to first order in 
IAVI/V 0 « 1. The suggestion here is to represent *(r) in a manner 
similar to p' (r) . Then 
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/(rk<r$r|<+i ( 0,0) = 0 (rk<rSrk+l)] jSj (0,^) = ^kjSj 
where ^kj is the jth spherical harmonic coefficient of the slowness in 
the kth layer for phase £(i). The path L(i) will depend upon jM(r), so 
the inverse problem is non-linear. I suggest both descending and 
ascending path segments of respective lengths L^ki and L a ki in layer k 
be determined initially from the mean state and subsequently from the 
results of the previous deep iteration. These ray path segment lengths 
are assigned to the mean location of the segment (flik.^k)- In this 
approximation, with the sum over k made explicit for descending and . 
ascending segments, 

= E [^kjSj (0kr^k)] Lki + T. tykjSj (0k>0k)] Lki 
k=l k=K 

i i i i 

- u-jhjSj (0u»^u) ~ vihjSj (0vi^v) 

where, to account for the effect of CMB topography on the path length, 
(i) uj = vj = 0 for paths not touching the CMB (e.g., P, S, PP, PS, 
etc.); (ii) ui = s 0 (r 0 = c + )/sint,u and vj = s 0 (r 0 =c + )/sini.v for 
reflections with incident angle o u and reflected angle at (^u.^u) = 
(0 v ,tv) (e.g., PcP, PcS, etc.); and (iii) ui = [s 0 (c _ ) - s 0 (c + )]/sini, u 
and vi = [s 0 (c _ ) - s 0 (c + ) ] / s i nt v for phases leaving the mantle with 
incident angle i, u at (0u,^u) and reentering the mantle with exit angle 
lm at (0 v ,0v) (e.g., PKP, PKS, etc.). Note and 6 V can be corrected 
using previous estimates of h. 

The foregoing expression for r i can be rewritten as 

d i a i u v 

r i = Z i j k^k j + Zijk^jk - Mijhj - Mijhj 

with the understanding that for seismic phase i on path L(i) between 
source r'(i) and receiver r'(i) Z d i j k (or Z a ijk) is the length of the 
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descending (or ascending) ray path segment in layer k times spherical 
harmonic j evaluated at the segment midpoint, while M u ij (or M v ij) is 
the travel time correction due to CMB topography at the point of core 
entry (or exit). This expression is more compactly represented by 
concatenating (i) the Zd and Z a tensors into Z, (ii) the matrices 
for P and S slownesses into and (iii) the M u and M v matrices into 
-M; then reorder the elements of f (and Z) into vector f (and matrix D) : 

r\ = DiijM + Mijhj . (12a) 

There remains the thorny problem of relating perturbation slowness 
to perturbation density. The differential slownesses for P-waves or S- 
waves are 

d^P = d(Vp)-l = -Vp _ 2dVp 

d*S = d(V S )-l = -V S -2dV S . 

With bulk and shear moduli K and p, V p2 = (K + 4p/3)/p and V$2 = p/p, so 

d p = V S -2d/4 - 2^V S -ldV S 

dp = Vp-2dK + (4/3)Vp-2d^ - 2/?Vp-ldVp 

which, upon linearization about the mean state (Vp 0 , Vso. Po> Ko, po) , 
are viewed as two equations in the three unknowns d p, dK, and dp. We 
can solve for 

dp = /jK-l[Vs 2 dK + 2/>VpVs(VsdVs - VpdVp)] 

where (Vp 2 - 4V$2/3) = K/p = (9p/9/>)ad is the adiabatic sound speed. 
Unfortunately, dp cannot be determined without additional information on 
dK (or dp) . 

The suggestion is to treat perturbation slowness as if directly 
proportional to perturbation density; however, the constant of 
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proportionality may vary radially (with the mean state). Lateral 
variations of this 'constant' are omitted for simplicity - as are 
anisotropies in the fourth rank-tensor of elastic constants and the 
complexities of attenuation. Then the perturbation density coefficient 
for spherical harmonic j in layer k is 

p -1 p s -1 s 

/>kj = [Ckk'l f^k ' j = [Ckk ' ] l^k'j 

or, in the vector notation, 

iM = cnvr • (12b) 

The elements of the diagonal matrix C are the different constants in 
each layer - be it mainly olivene; ol ivene-spinel; ferro-magnesian 
silicate perovskite; stishovite, non-stoiehiometric ferro-magnesio- 
wttstites, and ferro-silicides; or iron (Knittle & Jeanloz, 1989). 

If, within each layer, perturbations in temperature, pressure, and 
composition were directly proportional to density perturbations and if 
the partial derivatives of K with respect to temperature, pressure, and 
composition were laterally homogeneous, then (12b) would be fully 
justified. Hopefully, the information on dK needed for a more realistic 
equation of state can be obtained from mineral physics or mantle 
dynamics. Because d p, dp, and dK, hence dV$ and dVp, are caused by 
temperature, pressure, and composition perturbations associated with 
departures from the mean state of hydrostatic equilibrium with an 
adiabatic temperature gradient, both a perturbation equation of state 
and the equations of motion for the mantle will likely be needed. Then 
mantle circulation will have to be parameterized, the parameters 
included, and estimates thereof constrained to fit plate motion and 
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deformation data. (A parameterization of mantle circulation might also 
be tied to a parameterization of elastic anisotropy). 

In the interim, eqns. (12b) and (12a) allow the seismic portion of 
eq. (6) to be written 

[(Ti - r i) Wij (Tj - rj)] = [Ti - DiiCll'/>l' - Mj m h m ]Wij 

[Tj - DjlCllVr - Mj m h m ] 

= [T - DC/? - Mh]Tw S [T - DC p - Mh] (13) 

4. INVERSION 

With parameterization eqns. (9), (11), and (13), the constraint 
enforcing the geostrophic radial vorticity balance written (Bv)TA m (Bv) , 
and optional biases towards prior estimates of the fluid flow v°, 
topography h°, and density p° (possibly from 'disjoint' inversions), the 
objective function given by eq. (6) becomes 


* m * • er • 

A 2 = [g - Av]Tw [g - Av] + [D - QThv - q]Tw [Q - QThv - q] 

_ ge 

+ [c - F p - Hh]Tw [c - F p - Hh] 

+ [T - DC p - Mh]Tw S [T - DC p - Mh] 

+ (Bv)TAm(Bv) + (v - uO)Ta v (u - v°) 

+ (h - hO)TA h (h - ho) + (p - pO)fk p (p - po) . (14) 

In addition to the non-linearity explicit in the Earth rotation term 
(Qho = = Q v h) , recall that A depends upon 7, hence v; q depends upon 

w, hence both h and v; and D and M depend upon hence upon p and h. 

The minimization of A 2 is therefore a profoundly non-linear problem. 
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The attack on this problem offered here is based on iterative solution 
of the linearized problem. 

Initial estimates of A, q, D and M can be calculated from either 
the 'data' or the mean state. Initial estimates of Qh and Q w can be 
calculated from models of h and y obtained by 'disjoint' inversion. Use 
the initial estimates yi , hi, and pi with i = 1 (or 0 if necessary) to 
solve the forward problems of piecewise steady motional induction, 
changes in Earth rotation, geoid determination, and calculation of 
travel times. Such forward calculations give the predictions 7, w, f, 
and r. The differences between the 'data' g, 0, c, and T and these 
predictions are the residuals 5g, 50, 5c, and 5T. These residuals 
provide the 'data' for the first joint inversion (5g, 50, 5c, and 5T) . 
They also define the residual objective function, which is parameterized 
in terms of prior estimates ( y ° , h<\ ^0) ( initial values (yi , hi, pi), 
and parameter corrections (5yi , 5hi, dpi). The residual objective 
function (not shown) is minimal only if it is extreme, so parameter 
corrections can be estimated by setting the partial derivatives of it 
with respect to the parameter corrections equal to zero. 

In the linearized treatment, the elements of A, q, Qh, (] u , D, and M 
in the residual objective function are treated as if independent of the 
parameter corrections. Then the partial differentiation is easy and, 
to first order, the corrected parameters are the initial parameters 
(yi, hi, pi) plus 

— m - — er — - 1 

5v = {ATW A + QhTw Qh + BfA^B + Ay} 

{AlW^g + QhTw ef [5n-Q”5h-q] - BfA m Byi - A v [v 1 -v°] } (15a) 
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T er -r ge T S -1 

$h = {QvTw Qv + hTw H + MTw M + Ah} 



{QvTw [6Q-Qh6v-q] + (15b) 

HTw 9e (5c-F5/j) + MTw S [5T-DCty] - Ah[hi-ho]} 

ge s “1 

5p = {FTW F + [DC]TW [DC] + A s } 

T ge „ s 

{FTW [<5c-H<5h] + [DC]Tw [5T-M<5h] - A^-/^]} . (15c) 

Clearly, equations (15) are not fully reduced because of the terms 
Q v ffh (15a), Qhtfv and F 5p (15b), and H<5h and Mdh (15c) on the right. 
However, (15a) can be used to eliminate 6v from (15b) and (15c) can be 
used to eliminate 6p from (15b). The resulting expression can be solved 
for <5h. Then 5v and 6p can be determined by back substitution into 
(15a) and (15c). 

To make this clear, symbolically rewrite eqns. (15a-c) as 


6v = V-1[G - E5h] (16a) 
<5h = T-l [C - U5v - S 6p] (16b) 
6p_ = R-1[P - X5h] (16c) 
Then 

<5h = [T - UV-lE - SR-lX]-l[C - UV-lG - SR-lp] . (17) 


Substitution of eq. (17) into (16a) yields substitution of (17) into 
(16c) yields 6p. Clearly CMB topography is assigned the pivotal role. 

The corrected parameters can be used to solve the forward problem 
again and obtain new residuals. Then the corrections can be estimated 
again. After a satisfactory number of such "shallow" iterations, the 
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elements of A, q, Qv, Qh, D, and M can be recalculated using the 
predicted values 7(t), »(t) , £, and Vp(r) and Vs(r). Such recalculation 
is termed "deep" iteration; it might even include earthquake relocation. 
Then shallow iteration can be tried again, followed by another deep 
iteration, etc.. Shallow iteration is merely intended to provide a 
refined set of corrected parameters for use in the seemingly more 
burdensome recalculation of matrix elements required for deep iteration; 
it is considered optional for the small parameter corrections 
anticipated. Shallow iteration allows more use of residuals calculated 
by accurate numerical solution of the (time-dependent, non-linear) 
forward problem; such residuals ought not be replaced by coarse linear 
approximations. Clearly, the iteration process can be repeated until 
either an adequate fit is obtained or the Earth model is abandoned. 

Convergence of the iteration scheme is, of course, assured if the 
biases favoring prior estimates v°, ho, and p° (typically measured by 
by the diagonal elements of positive definite K v , Ah, and hp) are strong 
enough. Convergence is apparently neither prohibited nor guaranteed 
when either confidence in the prior estimates is eroded or when the 
prior estimates are replaced with the most recent estimate. In the 
latter case (which seems consistent with deep iteration), the (Ay.Ah.A^) 
might serve as convergence factors which keep the corrections so small 
as to avoid severe violation of the linearization when seeking small 
weighted residual variance. More sophisticated methods of non-linear 
optimization are possible, but lie outside the scope of this paper. 

Even in the linearized case, care is needed to avoid baseless bias, 
over-parameterization, and confusion of anticipated parameter error 
estimates (from the covariance) with the significance of the residuals. 
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5. DISCUSSION 

To obtain such a very simple deep Earth model complete through 
harmonic degree and order 10 with a 9-layer mantle, at least 1,440 
parameters need to be determined by iterative linearized least squares: 
120 coefficients representing the CMB topography, 240 coefficients 
representing the core velocity field per time interval, and 1,080 
coefficients describing the perturbation density. Symmetric storage of 
a 1440 x 1440 symmetric matrix requires but 10® words of computer 
storage. If more than one interval is considered, then more core flow 
coefficients are needed. For a 10*® degree, single-interval model with 
29 layers, there would be 3,840 parameters; symmetric storage of a 
3,840 x 3,840 matrix would require 7.4x10® words. Such matrices can be 
manipulated and, if wel 1 -conditioned, inverted with existing computers. 

The condition of the matrix depends upon the "data" selected, the 
assigned weights, the number of parameters estimated, and the confidence 
assigned to prior estimates. Note that truncation of the model need not 
suppose that higher degree parameters are zero; only that such unmodeled 
parameters contribute to the residuals and may contribute to model 
error. If the weights are to reflect covariance of unmodeled signal as 
well as data error covariance, expectations regarding unmodeled signal 
may be developed by studying the residuals obtained during numerical 
experimentation and, of course, the unmodeled processes themselves. 

In the very simple Earth model considered here, CMB topography 
provides the essential link between the diverse geophysical data types. 
Of course it is by no means clear that the connection between the 
magnetic side of the problem and either the geoid or the seismic sides 
of the problem is strong enough to warrant detailed calculations. This 
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connection is provided only through the three components of the decade 
fluctuations in earth rotation. Yet this connection can be strengthened 
by including the kinematic effect of CMB topography on core flow and 
thus the predicted secular geomagnetic variation. Furthermore, this 
connection might be strengthened upon parameterization of a more 
realistic Earth model, inclusion of more data types, and application of 
more physical constraints - as outlined in the introduction. 
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